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ABSTRACT 

Motivated by the interface model for the solar dynamo, this paper explores the complex 
magnetohydrodynamic interactions between convective flows and shear-driven insta- 
bilities. Initially, we consider the dynamics of a forced shear flow across a convectively- 
stable polytropic layer, in the presence of a vertical magnetic field. When the imposed 
magnetic field is weak, the dynamics are dominated by a shear flow (Kelvin-Helmholtz 
type) instability. For stronger fields, a magnetic buoyancy instability is preferred. If 
this stably stratified shear layer lies below a convectively unstable region, these two 
regions can interact. Once again, when the imposed field is very weak, the dynami- 
cal effects of the magnetic field arc negligible and the interactions between the shear 
layer and the convective layer are relatively minor. However, if the magnetic field is 
strong enough to favour magnetic buoyancy instabilities in the shear layer, extended 
magnetic flux concentrations form and rise into the convective layer. These magnetic 
structures have a highly disruptive effect upon the convective motions in the upper 
layer. 

Key words: convection - instabilities - (magnetohydrodynamics) MHD - Sun: in- 
terior - Sun: magnetic fields 



1 INTRODUCTION 

The 11 year solar magnetic cycle is driven by a hydro- 
magnetic dynamo. However, the exact nature of this dy- 
namo mechanism is still not fully understood, and there 
are several scenarios that seek to explain the obser ved be- 
havio ur. The well-known "interface" dynamo model (|Parkerl 
1 19931 ). is based on the idea that the dynamo operates 
in a region that straddles the base of the solar convec- 
tion zone and the stably stratified r egion that lies be- 
neath (for some recent reviews see lOssendriiven 120031: 
iProctoJ l200d : iDormv fc Sowardl 120071 : ISilverd 12008 ). Al- 
though this is a conceptually appealing model for the so- 
lar dynamo, the only numerical investigations of the in- 
terface dynamo have be en based upon mean-field dynamo 
theory (see, for example, Charbonncau & MacGrcgor 1997; 
Chan. Liao. Zhang fc Jones! I2004J : IZhang. Liao. SchubertJ 



been possible to demonstrate the operation of the interface 
dynamo by carrying out three-dimensional simulations of 
compressible magnetohydrodynamics. Given these compu- 
tational constraints, it makes sense to investigate different 
components of the interface dynamo in isolation. 

An important feature of the region below the solar 
convection zone is the solar tachocline (see Spiegel fc Zahnl 



2004 Bushbv 20061). In mean-field theory, several aspects 



of the dynamo model (particularly the effects of turbulent 
convection) are parametrised. However, the resulting coef- 
ficients are poorly determined by both theory and observa- 
tions. Due to the computational costs involved, it has not yet 
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Il992l ; Iciiristensen-Dalsgaard fc Thompson! [2007, and refer- 
ences therein), which takes the form of an intense radial 
gradient of the solar differential rotation. At the heart of 
the interface dynamo scenario is the idea that weak poloidal 
magnetic fields can be amplified by the intense shears in 
the tachocline, leading to the production of strong toroidal 
(azmimuthal) magnetic fields. In the standard interface dy- 
namo model, these poloidal magnetic fields are produced 
in the convection zone, and ar e pumped down into th e 
tachocline by the fluid motions (|Tobias et aZ.lll998l , l200ll ). 
In flux transport dynamo models, these poloidal fields are 
transported from the surface (to the tachocline region) b y 
a meridional circulation (see, e.g.. iDikpati fc Gilm an 2009). 
Wherever the poloidal field is generated, a mechanism is 
needed to produce flux structures that rise through the con- 
vection zone to the surface, where they emerge to form active 
regions. The most natural mechanism for inducing this verti- 
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cal tr ansport is magnetic buoyancy (|Parkerlll95l iNewcombl 
Il96j >. Strong coherent fields exert a magnetic pressure that 
leads to these magnetic regions becoming less dense than 
their surroundings. Provided that the ambient medium is 
not too stably stratified, instabilities can occur that would 
appear to allow strong fields to rise into the convection zone 
above. A full discussion of magnetic buoyancy and its im- 
portance in__relation to tachocline dynamics can be found in 
iHughed IHnl), while the role of the tachocline in the solar 
cycle is described bv lTobias fc Weiss! (|2007ft . 

Until recently, most studies have addressed the evolu- 
tion of magnetic buoyancy instabilities in a prescribed layer 



of m a gnetised fluid (see, for exampl e . ICattaneo fc Hughes 



19881 ; Matthews. Hughes fc Proctor! 19951; IWissink et al 



200ol : lFanll200ll ; iKersale. Hughes fc TobiasH2007l ). However, 
it is not immediately obvious that a realistic velocity shear 
can produce strong enough magnetic fields to become un- 
stable to buoyancy modes, particularly in the very stably- 
stratified tachocline. To become buoyant the fields must 
exert strong Lorentz forces, which will also retard the 
flow and resist the field amplification. The linear evolu- 
tion of magnetic buoyancy instabilities in a compressible 
magnetic layer, with an ali g ned v elocity shear, was con- 
sidered bv lTobias fc Hughes! (|2004h . They found that mag- 
netic buoyancy instabilities tended to be stabilised by a 
strong velocity shear. Recent numerical calculations have 
started to address the more complex problem of the non- 
linear evolution of shear -driven magnetic buoyancy instabil- 
ities (see, for example. [ Brummell. Cline fc Cattaneol [2002] ; 
Cattaneo. Brummell fc Clinell2006l ; IVasil fc Brummellll2008l . 
20091 ). 



Using a combination of high r esolution numerical simu- 
lation s and analytical calculations. I Vasil fc Brummell! ( 2008, 
2009) investigated the stability of a magnetic layer that 
is generated by the action of a strong vertical velocity 
shear upon an imposed uniform magnetic field. They ar- 
gued that no magnetic buoyancy instability would be pos- 
sible, in the stably-stratified tachocline, unless the mag- 
nitude of the velocity shear were many orders of mag- 
nitude larger than the inferred radial shear. This would 
have profound consequences for the hydrodynamic stabil- 
ity of the shear. Defining the Richardson number, Ri, to be 
the square of the Brunt- Vaisala frequency divided by the 
square of the velocity gradient, a necessary condition for 
hyd rodynamic stability i s that Ri > 1/4 (see, for exam- 
ple, [Chandrasckhar 1961). In the tachocline, the Richard- 
son number is estimated to be many orders of magnitude 
larger than that given by this stability bound, which im- 
plies that the shear is stable. However, if the velocity shear 
is strong enough that the stabi lity condition is n o t sat- 
isfied , as in the calculations of IVasil fc Brummell l|2008l . 
2009), then this system will be subject to shear instabil- 
ities (of "Kelvin-Helmholtz" type). Clearly the situation 
becomes more complicated in the presence of an imposed 
magnetic field, and the subsequent evolution depends cru- 
cially (and highly non-trivially) upon the strength of this 
magn etic field. This is an int eresting problem in its own 
right. iHughes fc Tobias] (120011 ) considered the linear evolu- 
tion of magnetised shear instabilities, whilst the nonlinear 
proble m has also been studied in unstratified compressible 
fluids (iFrank. Jones. Ryu. Gaa laas 199£ ;lRvu. Jones. Frank] 



|2000| ; iPalotti, Heitsch. Zweibel. Huand I2OO8I ) as well as 



isothermal stratified layers jBriiggen fc Hillebrandtl l200ll ). 
The recent review article by iGilman fc Callvl \200H) de- 
scribes global magnetohydrodynamic shear instabilities in 
the tachocline. 

Although the hyd rodynamic stabilit y of the vel ocity 
shear was discussed by I Vasil fc Brummelll (120081 , |2009j) , the 
most important idea in their work was the suggestion that 
shear-driven magnetic buoyancy instabilities can only oc- 
cur at very small values of the Richardson number. Since 
this would appear to be incompatible with the tachocline, 
this would have dire consequences for solar dynamo mod- 
els. However, these results are not conclusive. Firstly, their 
calculations were all performed using a fixed value for the 
imposed magnetic field strength. This is clearly an im- 
portant parameter, since the Lorentz force plays a cru- 
cial dyjiamjcalrol^More importantly, recent calculations 
by ISilvers et al\ (120091 ) have confirmed, as already known 
for th e onset of magnetic buoyancy without shear (|Hughed 
120071 ). that the onset of magnetic buoyancy instabilities de- 
pends upon the ratio of the mag netic to thermal di ffusiv- 
ities. At high Reynolds numbers, ISilvers et al\ ((2009) have 
shown that magnetic buoyancy instabilities can be excited 
with a weaker (hydrodynamically-stable) shear if the ther- 
mal diffusivity is much greater than the magnetic diffusivity 
(s omething that was not th e case in the original calculations 
of lVasil fc Brummelll l|2008h ). This is a more encouraging re- 
sult from the point of view of the solar dynamo, although 
more work remains to be done. 

Clearly, the parametric dependence of the instabilities 
of a forced shear flow, in the presence of a magnetic field, 
is still not fully understood. One of the aims of this pa- 
per is to enhance our understanding of these instabilities 
via a partial exploration of parameter space. In particular, 
we focus attention upon the effects of varying the strength 
of the imposed magnetic field (although some variations in 
other parameters are also considered). Once the evolution of 
this system has been studied in a single convectively-stable 
layer of compressible fluid, we move on to consider a more 
complicated "composite" model, which combines the stably- 
stratified shear layer with an overlying convectively-unstable 
region. This composite model enables us to address the in- 
teresting question, so far largely unexplored, of how buoy- 
ant magnetic flux might interact with the fluid in the lower 
convection zone. In order to limit the computational ex- 
pense of this parametric survey, we choose a stronger veloc- 
ity shear than that considered bv lSilvers et al. (2009). This 
enables us to drive buoyancy instabilities at lower Reynolds 
numbers, which means that fully resolved numerical simula- 
tions can be carried out with a coarser numerical grid. Al- 
though our chosen flow is hydrodynamically unstable, it is 
sti ll much weaker than the target shear that was considered 
by I Vasil fc Brummelll (120081 . 120091 ). being mildly subsonic as 
opposed to highly supersonic (though still much stronger 
than that found within the tachocline). Although we are 
not exploring the rather extreme parameter regime that is 
directly relevant for the tachocline, our choice of parameters 
allows us to enhance our basic understanding of the interac- 
tions between magnetic buoyancy and convective instabili- 
ties. Future research (which will rely upon this work) will 
focus upon these phenomena at higher Richardson numbers. 

The plan of the paper is as follows: In the next section 
we describe the set up of the model problem. In section 3, we 
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present numerical results from this model, describing the in- 
teractions between magnetic and hydrodynamic instabilities 
in a single stably-stratified polytropic layer. In the follow- 
ing section, we describe the (more complicated) problem of 
shear-driven instabilities in the composite model. Finally, in 
section 5, we conclude with a discussion of the astrophysical 
significance of these results. 



2 THE MODEL 

The model we use is similar to that of several pre- 
vio us studies of magnetoconvection (see, for exam- 
ple, IMatthews. Proctor fc Weissl ll 995l ; iBushbv fc Houghton! 
l2005l ; iLin. Silvers fc Proctonl2008h . We consider the evolu- 
tion of a plane layer of electrically-conducting fluid, which 
is heated from below, in the presence of a magnetic field 
that is initially uniform and vertical. Accordingly, we adopt 
a Cartesian frame of reference such that the z-axis points 
vertically downwards, parallel to the constant gravitational 
acceleration, g. Defining d to be some characteristic length- 
scale (e.g. the depth of the convection zone in the com- 
posite model), this fluid occupies the region ^ x,y ?C 
Ad and ^ z ^ dzo. We set zo = 1 in the single 
layer calculations that are described in the next section 
(in which case d corresponds to the layer depth), whilst 
zo = 3 in the composite model. Varying the parameter A 
enables us to change the width of the domain without al- 
tering the value of d. In contrast to most previous stud- 
ies of_jmyinetoc^^ the recent paper 
bv lKapyla. Korpi fc Brandenburg! ; 2008). we investigate the 
evolution of this system under the influence of a forced hor- 
izontal shear flow in the x-direction. 

Throughout this paper, we assume that the shear vis- 
cosity, n, the magnetic diffusivity, rj, the permeability of free 
space, /io, and the specific heat capacities at constant pres- 
sure and density (cp and cv respectively) are all constant 
properties of the fluid. The thermal conductivity, K(z), is 
assumed to be a function of z. Defining p to be the fluid 
density, T to be the temperature, u to be the fluid velocity 
and B to be the magnetic field, the governing equations for 
the evolution of this compressible fluid are given by: 



dui 
dxi 



+ 



dxi 



2 duk 

3 dx k 



(7) 



| + V. ( pu) = 



(1) 



^ + (u-v)u - -vr + mz-fV- ( -,,<-. ix] (■>) 



dT 



dB 
dt 
V B 



+ — (V x B) x B + V- ((mt) 

= -PV-u + V-[if(z)VT] (3) 

77IV x B| 2 ur 2 
+ L + V 

= V x [u x B — 77V x B] (4) 
= 0, (5) 



where the pressure P satisfies the perfect gas law 

P = R* P T, (6) 

(defining _R» to be the gas constant) and the components of 
the viscous stress tensor t satisfy 



Finally the scalar quantity, Uo(z), represents the horizontal 
shear flow. The corresponding forcing term in Equation [2] 
ensures that any imposed shear of this form is a solution 
of the horizontal component of the momentum equation (in 
the absence of any other motions). 

The boundary conditions for these variables are consis- 
tent with those of an idealised model. All variables are as- 
sumed to satisfy periodic boundary conditions in the x and y 
directions. The upper and lower bounding surfaces (at z = 
and z = dzo respectively), are assumed to be impermeable 
and stress-free, and it also assumed that the magnetic field 
is vertical at these boundaries. The upper boundary is held 
at fixed temperature, whilst the heat flux passing through 
the lower surface is assumed to be constant. This implies 
that: 



du x 
dz 
du x 
dz 



du y 
dz 



dz 



- B x 
B x 



B y = 0, T = T at z = 
dT 



B v =0 



dz 



(8) 

C at z = dz , (9) 



where C is a constant that will depend upon the initial con- 
ditions of the model. Note that the choice of a stress-free 
boundary condition for u x implies that the imposed shear, 
dUo(z)/dz, should also be zero at these surfaces. 

These equations can be expressed in non- 
dimensional form, using the s calings described by 
IMatthews. Proctor fc Weiss! l| 19951 ). Lengths are scaled 
by d, whilst the density and temperature are scaled by their 
initial values at the top of the layer (pa and To respectively). 
Velocities are scaled in terms of the isothermal sound speed 
at the upper surface, ^/R,Tq, which suggests a natural 
scaling for time of d/^/RTTo- Magnetic fields are scaled in 
terms of the strength of the initial vertical magnetic field 
Bo- Finally, we define Ko to be the value of K(z) at the 
upper surface. When these scalings are substituted into 
the governing equations, we obtain several non-dimensional 
par ameters that are essentially i dentic al to those described 
by IMatthews. Proctor fc Weiss! l|l995l ). These include the 
dimensionless thermal diffusivity, re = Ko/dpoCp^J (R,To), 
the ratio of specific heats, 7 = cp/cy, the Prandtl number, 
a — pcp/Ko, and the ratio of the magnetic to the thermal 
diffusivity at the top of the layer, fo = rjCppo/Ko- Finally, 
F — Bq J RtTopoPo is the ratio of the squared Alfven speed 
to the square of the isothermal sound speed at the top 
of the layer. This parameter determines the dynamical 
influence of any imposed magnetic field. 

With appropriate choices for these non-dimensional pa- 
rameters, we solve the equations numerically using a parallel 
hybrid finite-difference/pseudo-spectral code. In this code, 
time-stepping is carried out with an explicit third-order 
Adams-Bashforth scheme. Horizontal derivatives are eval- 
uated in Fourier space, whilst vertical derivatives are calcu- 
lated using fourth-order finite differences (upwinded deriva- 
tives being used for the advective terms). In order to carry 
out these simulations, grid resolutions of 128 x 128 x 200 
mesh points were used for the single layer cases, whilst 
256 x 256 x 300 mesh points were used for the composite 
model. Some calculations were also carried out at lower spa- 
tial resolution and comparisons of the different resolutions 
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Param. Description Value 



a 


Prandtl Number 


0.05 


ni 


Polytropic Index 


1.6 


e 


Thermal Stratification 


0.5 


7 


Ratio of Specific Heats 


5/3 


K 


Thermal Diffusivity 


0.01 


F 


Magnetic Field Strength 


Variable 


Co 


Magnetic Diffusivity 


0.2 



Table 1. Parameter values for the single layer calculations. 



show that the instabilities and structures that emerge are 
physical and not artefacts of the discretization. 



3 THE SINGLE LAYER 

In this section, we consider the evolution of this system with 
the simplest possible initial configuration. An understand- 
ing of this system will help us to interpret the results from 
the next section, which deals with a much more complicated 
model problem. Throughout this section, we define the com- 
putational domain by setting A = 4 and zo = 1- Recalling 
that all lengths are scaled in terms of d, this implies that 
< x, y < 4 and < z < 1. 

3.1 Parameters 

The behaviour of this model depends crucially upon the ini- 
tial conditions that are imposed. The simplest non-trivial 
case to choose is that of a polytropic layer with a con- 
stant thermal conductivity, K(z) = Ko- In this case, the 
initial conditions are completely determined by two non- 
dimensional parameters, namely the (dimensionless) tem- 
perature difference between the upper and lower boundaries, 
9, and the polytropic index, m = gd/ 'R*Tq9. Neglecting the 
effects of viscous heating, it is straightforward to show that 
the governing equations have the following (dimensionless) 
equilibrium solution: T = (1 + 6z), p = (1 + 9z) m , B z = 1, 
u x = Uo(z), Uy = u z = B x = By = 0. Of course the effects 
of viscous heating will lead to a departure from this equilib- 
rium, but we have verified (by direct calculation) that the 
departure is negligible over the time-scales that are consid- 
ered in this paper. Therefore the above "equilibrium" solu- 
tion (together with a small, random, thermal perturbation) 
is used as an initial condition for all the simulations that are 
described in this section. Note that these initial conditions 
imply that the lower boundary condition for temperature 
(see Equation O becomes dT/dz = 9 at z — 1. 

This system of equations has a large number of di- 
mensionless parameters, making it impractical to conduct 
a complete survey of parameter space. Therefore we focus 
primarily upon varying the strength of the magnetic field, 
holding all other parameters fixed (although a small num- 
ber of runs with different parameter values were also car- 
ried out). The parameter choices are summarised in Ta- 
ble [1] Note that this choice of 9 implies that the layer 
is weakly stratified. Setting m — 1.6 and 7 = 5/3 im- 
plies that stratification in the layer is mildly subadia- 
batic. This choice of parameters is appropriate for the 
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Figure 1. Top: u x as a function of depth for the single layer 
model. Bottom: du x /dz as a function of depth for the single layer 
model. 

stably stratified solar tachocline (IVasil fc Brummelll 120081 : 
IChristensen-Dalsgaard fc Thompsonl I2007T ). Note also that 
the parameter values that are given in Table [T] imply that 
the dynamical effects of the magnetic diffusivity, the viscos- 
ity and the thermal diffusivity are much more significant in 
this model than in the solar interior. This is because the 
dissipative length scales associated with the solar interior 
can not be resolved using any current computer. However, 
by setting 1 > C,o > u, we ensure that these are ordered in 
the same way as in the solar interior, i.e. the thermal dis- 
sipative cutoff scale is larger than the magnetic dissipative 
cutoff scale, which is in turn taken to be larger than the 
viscous scale. 

Finally, we must specify a suitable initial shear flow for 
this system. We set 

U {z) = 0.577(1 + tanh[20(z - 0.5)]) (10) 

as shown in Figure QJtop). The hyperbolic tangent gives a 
smooth velocity field, varying from u x — at z — up to 
u x = 1.154 at z = 1. The width of the shear region is suffi- 
ciently small that the departure from a stress-free condition 
at the boundaries is comp arable with the numerica l error of 
the scheme. The results o f IVasil fc Brummelll (| 20081 ) suggest 
that a stronger shear promotes magnetic buoyancy instabil- 
ities. We have maximised the shear velocity subject to the 
constraint that the horizontal flow speed never exceeds the 
adiabatic sound speed, whilst also ensuring that the peak 
mach number of the flow is identical to the peak mach num- 
ber of the shear in the composite model (see the next sec- 
tion). Note that the fluid Reynolds number of this shear 
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(based upon the peak velocity and the width of the shear 
layer) is appro ximately 300, which is much smaller than in 
other studies (|Vasil fc Brummelll l2008t 120091 : ISilvers et all 
2009). As discussed in the Introduction, this enables us 
to carry out fully resolved simulations with comparatively 
modest numerical grids. 

3.2 Results 

Having set up the model problem for a stably-stratified poly- 
tropic layer, we now investigate the effects of varying the im- 
posed magnetic field. In order to achieve this, we carry out 
a series of numerical simulations for different values of the 
parameter, F (keeping all other parameters fixed, as shown 
in Table [TJ. In the absence of a magnetic field (F = 0), 
we find that the system is unstable to a shear flow (Kelvin- 
Helmholtz type) instability. Initially, we see rolls forming in 
the x — z plane, as shown in Figure [2] This then rapidly 
evolves into a three-dimensional time-dependent flow. The 
implications of imposing a magnetic field across this compu- 
tational domain depend upon the strength of the imposed 
magnetic field. If the initial field is weak (say F — 1/90000), 
then the effect of the field on the evolution of the instability 
is negligible, as the solutions are virtually indistinguishable 
from the hydrodynamic case. The magnetic field is simply 
advected with the resultant flow as a passive vector field. 
However, if we increase the strength of the magnetic field, 
we find that it starts to have a dynamical influence. Fig- 
ure [3] shows the evolution of a simulation with a magnetic 
field strength determined by F = 1/9000. In this case, the 
magnetic field does reduce the vigour of the instability, lead- 
ing to more ordered motions (particularly at early times). 
This behaviour is easy to explain. The shear flow instability 
acts so as to bend the magnetic field lines parallel to the 
velocity shear. A strong field tends to resist this process, 
thus inhibiting the instability. However even in this simula- 
tion, as in the hydrodynamic case, the instability eventually 
develops three-dimensional structure. 

The character of the instability changes further as we 
increase the strength of the imposed magnetic field. Results 
for F — 1/900 are shown in Figure [4] The resulting hor- 
izontal magnetic field is now strong enough to completely 
suppress the shear flow instability. Rather than generating 
fluctuations parallel to the shear, the initial instability is 
a short wavelength interchange instability, with almost all 
variation (at least initially) in the y direction. These inter- 
change modes are typical of a m agnetic buoyancy instability 
l|Newcomblll96ll ; [Hughes 1985:). During these early stages, 
the developing structures are similar to those found in two- 
dimensional calculations of t he break up of a magnetic layer 
in the absence of a shear l|Cattaneo fc Hu ghes 1988]). At 
later times some longer wavelength variation in the x direc- 
tion does appear - this three-dimension al evolution is similar 
to that found bv lWissink et al\ (|200Ch . 

While the focus of this section is to explore the effects of 
varying the magnetic field strength, we note that there are 
other parameters that can be varied (subject to computa- 
tional constraints). A reduction in the Prandtl number leads 
to a reduction in the viscous dissipation relative to the other 
diffusivities. This also increases the fluid Reynolds number. 
We carried out runs with lower values of the Prandtl num- 
ber, and found that reducing this parameter by up to a fac- 






Figure 2. Density perturbation snapshot for the hydrodynamic 
case (F = 0) at times t = 2.59 (top), i = 5.20 (middle) and 
t = 7.69 (bottom) 
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Figure 3. Density perturbation snapshot for F=l/9000, at 
t=5.21 and 7.78. 



tor of ten (going down to a ~ 0.005) appears to have little 
effect upon the vigour of the instability at F = 1/900. This 
suggests that there is very little dependence upon the fluid 
Reynolds number in this parameter regime. The effects of 
varyi ng both k and a ha ve not been systematically studied 
here; ISilvers et all (2009) have shown that k in particular 
can play an important role for larger values of the Richard- 
son number, and investigation of a full range of diffusivity 
ratios will be the subject of future work. 



4 THE COMPOSITE MODEL 

In this section, we consider a more complicated model prob- 
lem, consisting of a piecewise polytropic atmosphere. This is 
intended to be a highly idealised representation of the region 
straddling the base of the solar convection zone. In order to 
achieve this, we consider a deeper computational domain, 
corresponding to zo = 3. We also choose a wider computa- 
tional domain, by setting A = 8. Recalling the definitions of 
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0.05 




Figure 4. Density perturbation snapshot for F=l/900, at 
t=10.93 and 13.88. 



these parameters, this implies that the computational do- 
main is defined by ^ x, y < 8 and < z < 3. 



4.1 Parameters 

Other than the dimensions of the computational domain, 
the main difference between these calculations and those of 
the preceding section is that the polytropic index of the do- 
main is now a function of depth. We split up the domain into 
three layers of unit depth. In the top layer (0 ^ z ^ 1), we 
choose a polytropic index of mo — 1, which implies that this 
region is convectively unstable. Like the single layer from 
the previous section, the middle region (1 ^ z ^ 2) is con- 
vectively stable with mi = 1.6. The lower layer (2 ^ z ^ 3) 
is also convectively stable but with a much larger polytropic 
index, ni2 = 4. The primary purpose of the lower layer is to 
lessen the impact of the rigid lower boundary. Any descend- 
ing convective plumes that reach z = 2 can simply pass into 
the lower region without "splashing" back and interfering 
with the other dynamics in the system. In order to achieve 
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Figure 5. Top: u x as a function of depth for the composite layer 
model. Bottom: du x /dz as a function of depth for the composite 
layer model. 



Figure 6. The initial temperature and density profiles for the 
composite model. 



Par am. 


Description 


Value 


(7 


Prandtl Number 


0.05 


mo, mi , ni2 


Polytropic Indices 


1.0, 1.6, 4.0 


7 


Ratio of Specific Heats 


5/3 


K 


Thermal Diffusivity 


0.0385 


F 


Magnetic Field Strength 


variable 


Co 


Magnetic Diffusivity 


0.1 



Table 2. Fixed Parameter Values 



the required piecewise-polytropic structure, we choose the 
following depth-dependent thermal conductivity profile: 



K(z) 



tanh 



+ 



Ko 
2 

Kp (m 2 + 1) 
2 (mo + 1) 

Kg (mi + 1) 
2 (m + 1) 



V 0.02 )_ 
1 + tanh 



1 — tanh 



(11) 



2-2 

0.02 



{ -\ tanh ( 

\ 0.02 J \ 



z - 1 
0.02 



where Ko = K(0) (as before). The tanh profiles ensure that 
the conductivity varies smoothly between each region. 

In this composite model, our aim is to investigate the ef- 
fects that any shear instabilities in the mid-layer have upon 
an established pattern of convection. Therefore, we only in- 
troduce the shear once the convection in the upper layer 
has become fully developed. This is achieved by integrating 
the equations without any horizontal forcing until t « 40. 
The shear is then introduced at this point, along with the 
corresponding forcing term in Equation [2] Once it has been 
introduced, the shear has the same structure as that in the 
single layer model but is now centred at the mid-plane of 
the middle region (at z = 1.5). Thus the imposed shear now 
has the form: 



U {z) = 1 + tanh[20( 2 - 1.5)], 



(12) 



as shown in figure [5] Note that the amplitude of the shear 
is chosen so that the local Mach number of the flow is the 
same as for the single-layer case. 

In the absence of any imposed shear at t = 0, the ini- 



tial conditions for this model differ slightly from those in 
the single layer case. Here we choose a magnetohydrostatic 
initial condition, setting B z = 1 and u x = u y — u z = 
B x — By — 0. The equilibrium profiles for p(z) and T(z) 
are found numerically, and are shown in Figure [6] Note that 
the choice of the thermal boundary condition at the lower 
surface determines the extent of the thermal stratification. 
Setting dT/dz = 0.8 at z — 3 ensures that the temperature 
increases by 50% across the middle layer, as was the case for 
the single layer. 

The parameters for this composite model are chosen 
so that the conditions in the middle layer are as similar as 
possible to those for the single layer calculation. Note that 
this requires some rescaling of k and (o- These parameters 
are shown in Tabled 



4.2 Results 

As for the single layer calculations in the previous section, 
we carry out a series of numerical simulations for different 
values of the imposed magnetic field (as measured by F). In 
addition to any effects upon the hydrodynamic instabilities 
of the shear, increasing F also reduces the vigour of any 
convective motions in the upper region of the computational 
domainQ The range of F is carefully chosen so that we cover 



1 Note that estimates from linear theory suggest that a value of 
F of approximately 0.5 is needed in order to completely suppress 
convective flows in the upper layer. Therefore, we are not near 
the convective stability boundary. If we were to increase F fur- 
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0.8 




Figure 7. Three-dimensional plots of the F = 0.0001 case in the 
absence of any velocity shear. Top: The temperature perturbation 
on the side of the computational domain and close to the upper 
surface at t m 40. Bottom: The same plot at t fs 51. 

the same values for the mid-layer plasma beta (the ratio of 
the gas pressure to the magnetic pressure) that were covered 
in the single-layer case. 

Initially, we set F — 0.0001, which corresponds to a 

ther, we would expect to see a transition to an oscillatory mode 
of convection before we reach the regime in which convection is 
completely inhibited. 




Figure 8. Three-dimensional plots of the F = 0.0001 case at 
t m 51, after the shear was introduced at t « 40. Top: The tem- 
perature perturbation on the side of the computational domain 
and close to the upper surface. Bottom: The horizontal magnetic 
field on the same surfaces. Note that in this figure, B x is nor- 
malised with respect to the imposed field, oc y/~F. 

weak imposed magnetic field. Given the relative complexity 
of this system, we first explore the dynamics that occur in 
the absence of any velocity shear. This case is illustrated in 
Figure [7] which shows the resulting pattern of convection 
at t « 40 (top) and t « 51 (bottom). In the convectively 
unstable upper layer there is a time-dependent cellular con- 
vective pattern consisting of warm, broad upfiows (which 
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correspond to the brighter regions in figure [7} surrounded 
by a network of narrower (darker) downflows. There is some 
modest convective overshooting into the middle layer. Fig- 
ure [8] also shows a snapshot of this system at (t ~ 51) 
but, in this case, the shear is introduced at t « 40. As 
was found in the single-layer case, the shear is subject to 
a Kelvin-Helmholtz type shear instability. This rapidly de- 
velops three-dimensional structure, producing perturbations 
that spread throughout the stably-stratified domain, pene- 
trating into the convective layer. Comparing Figure [7] and [8] 
we see that there is little evidence of anisotropy in the con- 
verting region. However we note that there is a significant 
influence on the convective transport of heat in the top part 
of the box, which is apparently due to the influence of the 
shear instability on heat transport in the middle layer. 

As the field strength is increased from this (effectively) 
kinematic level, the solutions follow a similar trend to the 
single-layer case. When F = 0.001, the dynamics in the mid- 
layer are similar in form to those shown in Figure|3] (although 
the overshooting convection from the upper layer adds some 
additional complexity to the resulting flows). Therefore, the 
dominant instability is still of Kelvin-Helmholtz type rather 
than a magnetic buoyancy instability, although magnetic ef- 
fects are starting to play a dynamical role. Interestingly, as 
the shear-driven motions from the stable layer interact with 
the convective layer, there appears to be a slight tendency 
for an elongation of the convective cells in the direction of 
the shear. This is a phenomenon that becomes more pro- 
nounced as F is increased. 

Increasing the field strength still further, so that F = 
0.01, we find that the dynamics change dramatically. This 
case is illustrated in Figures l9l and [lOl The imposed vertical 
magnetic field is now strong enough to reduce the vigour 
of the convection, though it has little effect on the hori- 
zontal scales of motion. Once the shear is introduced, the 
evolution is dominated by the shear-driven instabilities at 
the mid-layer. As in the single layer case, a transition has 
occurred so that the dominant instability is now magnetic 
buoyancy. Initially, this buoyancy instability takes the form 
of a two-dimensional (interchange) mode, although it soon 
develops three-dimensional structure, forming arching re- 
gions of magnetic flux that rise up through the convective 
upper layer of the domain. As these magnetic regions reach 
the upper layers, we see some concentration of the verti- 
cal magnetic flux, which forms localised concentrations near 
the horizontal boundaries of these rising features. The subse- 
quent motion is now strongly anisotropic, producing convec- 
tive cells that are predominantly aligned with the direction 
of the shear and the buoyant horizontal magnetic flux con- 
centrations. We note that the introduction of a shear flow at 
t ss 40 again leads to larger temperature deviations in the 
convectively-unstable region at later times (compared with 
the unsheared case). However this effect seems to become 
less pronounced as the field strength is increased. We at- 
tribute this phenomenon to the fact that there is less mixing 
in the stronger field regime, where there are larger structures 
present than in the weaker field cases. 

The phenomena discussed above can be related to pre- 
vious work on isolated buoyant flux tubes rising throug h 
the convection zone (see, for example. I Jouve fc Brunll2007h . 
Even though the tubes are isolated they have been shown to 
interact strongly with the convective flow. The present prob- 



9.6 




Figure 9. Three-dimensional plots of the F = 0.01 case in the 
absence of any velocity shear. Top: The temperature perturbation 
on the side of the computational domain and close to the upper 
surface at t S3 40. Bottom: The same plot at t ~ 51. 

lem is different in that the magnetic field is initially vertical. 
However, the shear creates a strong horizontal magnetic field 
and so the ultimate configuration is not dissimilar. 

5 CONCLUSIONS 

In this paper, we have presented some novel calculations to 
investigate the ways in which an imposed magnetic field in- 
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Figure 10. Three-dimensional plots of the F = 0.01 case at 
t as 51, after the shear was introduced at t « 40. Top: The tem- 
perature perturbation on the side of the computational domain 
and close to the upper surface. Bottom: The horizontal magnetic 
field on the same surfaces. Note that in this figure, B x is nor- 
malised with respect to the imposed field, <x yF. 

teracts with a shear flow in a convectively-stable layer (both 
with and without an overlying convective region). This in- 
vestigation was motivated by conditions at the base of the 
convection zone, and is relevant to the interface scenario 
for the solar dynamo. The most important interactions be- 
tween the shear layer and the convective region are due to 
the rising plumes that are induced by magnetic buoyancy. In 



the solar context the tachocline region, where the shear is 
expected to reside, is very stably stratified and there are 
questions regarding the efficacy of magnetic buoyancy i n 
this situation (see, for example. IVasil fc Brummelll [2009). 
Nonetheless we know that buoyancy instabilities do occur 
in the Sun, and so it seems worthwhile trying to understand 
aspects of their evolution, even though the correct parame- 
ter range might not yet have been reached. In this context, 
it is also worth noting that a diffusive instability, which is 
effective only when the thermal diffusivity is much higher 
than in the present paper, appears to allow for buoyancy- 
induced motion even when the shear is hydrodynamically 
stable, according to the Richardson number criterion. This 
mechanism been discusse d for several decades (see, fo r ex- 
ample, the discussions in lGilman]|l970l ; iHug hes 200?]) but 
has only recently been dem onstrated numerically in our ge- 
ometry IjSilvers el aZ.ll2009f ) . 

In our "strong-shear" parameter regime, there is an in- 
stability of the shear (of Kelvin-Helmholtz type) even when 
there is no magnetic field. This instability leads initially to 
perturbations that are primarily in the x — z plane. This in- 
stability subsequently develops structure in the y direction. 
We find that, for a sufficiently strong magnetic field, the hy- 
drodynamic instability is suppressed, leading to a magnetic 
buoyancy instability with strong variations in the y direc- 
tion, i.e. the direction perpendicular to both gravity and the 
shear flow. For calculations of the "composite model", the 
most important interactions between the shear layer and the 
convective region are due to the rising plumes that are in- 
duced by this magnetic buoyancy instability. This instabil- 
ity generates strong horizontal concentrations of magnetic 
flux that then rise through the convective layer. For the 
strongest fields surveyed, the dynamical effects of these flux 
concentrations are so significant that they destroy the con- 
vective pattern that would normally exist in the absence of 
a m agnetic field. As wit h isolated flux tubes (see, for exam- 
ple, |jouw^^run|[20o3) they effectively push the convective 
motion aside. It is important to note that even our strongest 
imposed field has a relatively weak effect upon the horizontal 
scales of convection that occur in the absence of shear. 

Our results are encouraging in that they show that 
there exists the possibility of inducing strong buoyant mag- 
netic flux structures through the action of horizontal shear. 
However they are preliminary calculations in the sense that 
they do not allow the buoyant flux structures to rise very 
far. We intend to perform further calculations, with a much 
deeper convective zone, so as to understand the later evo- 
lution of the rising flux structures. We also note that there 
is a swift transition as we vary our parameters between es- 
sentially passive structures and ones that strongly disrupt 
the convection layer. We intend to carry out a more exten- 
sive investigation of intermediate parameter ranges in which 
the role of downward pumping is import ant in counteracting 
the buoyancy effects (jTobias et aZ.|[200ll 'l . Finally, we intend 
to explore the interactions between shear-driven buoyancy 
instabilities and convective flows at higher Richardson num- 
bers (with a hydrodynamically stable shear). This will be a 
challenging problem to tackle numerically (requiring high 
numerical resolution), but results from these preliminary 
calculations constitute a firm foundation for future work. 
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